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^ Abstract 

O The Lyapunov exponents of a chaotic system quantify the exponential divergence 

of initially nearby trajectories. For Hamiltonian systems the exponents are related 
to the eigenvalues of a symplectic matrix. We make use of this fact to develop a new 
method for the calculation of Lyapunov exponents of such systems. Our approach 
avoids the renormalization and reorthogonalization of usual techniques. It is also 
easily extendible to damped systems. We apply our method to two examples of 
physical interest: a model system that describes the beam halo in charged particle 
beams and the driven van der Pol oscillator. 
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Chaotic dynamics has been observed in a wide variety of systems, including bio- 
logical systems, weather models, mechanical devices, plasmas, and fluids, to name a 
few. In the chaotic regime these systems exhibit exponential divergence of initially 
nearby trajectories, as embodied by the Lyapunov exponents of the system [I]. Over 
the past two decades or so there has been "intense activity" [2] directed toward the 
computation of these exponents resulting in several different numerical approaches 
[I] [3] [4]. The two difficulties associated with the computation of Lyapunov exponents 
are: (I) exponential growth of the separation vector (between the fiducial and nearby 
trajectories) and (2) the exponential collapse of initially orthogonal separation vectors 
onto the direction of maximal growth. Most conventional methods overcome these 
hurdles by intermittent numerical rescaling and reorthogonalization (through, e.g., 
the Gramm-Schmidt procedure [3]). Many chaotic systems are Hamiltonian or they 
can be transformed into a Hamiltonian system by suitable manipulations. However, 
none of the above general methods are designed to take advantage of this fact. 

The dynamics of classical Hamiltonian systems has an underlying symplectic struc- 
ture [5]. In recent years symplectic methods have been applied with great success to 
classical dynamical problems. The field of accelerator dynamics has been revolution- 
ized by the introduction of nonlinear symplectic maps as represented by Lie transfor- 
mations [6] [7]. Very long time integration of charged particle and planetary systems 
has been aided by the development of high order symplectic integration algorithms 
[8]. In this Letter we exhibit a symplectic map-based approach towards the calcu- 
lation of Lyapunov exponents. As shown below, this approach obviates analytically 
the need for rescaling and reorthogonalization in the numerical computation of the 
exponents. We have successfully applied our method to several systems, including the 
Duffing oscillator, the damped driven pendulum, the driven double-well system, the 
beam halo system, and the van der Pol oscillator (details will be presented elsewhere 
[9]). In this Letter, we will briefiy describe our general approach and expose in some 
detail the latter two examples. 

Consider a 2-m dimensional continuous-time dynamical system governed by the 
equations 



where z = {zi, ■ ■ ■ , Z2m) and similarly for F. Let Zq denote some given fiducial 
trajectory, and suppose we wish to study nearby trajectories. To do this, define 
deviations from the fiducial trajectory by letting Z = z — Zq, and linearize the above 
equations. The new set of equations for the deviation variables is 
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Our approach can be used whenever this hnearized set of equations can be derived 
from a Hamiltonian. From now on we will suppose that this is the case, and that 
we can write Z = {qi,q2, - ■ ■ ,qm,Pi,P2^' ' ' ^Pm)^ where qi and pi denote canonically 
conjugate coordinates and momenta, respectively. It follows that 

f^-WZ), (3) 

where {,} denotes the Poisson bracket, and where H is a homogeneous quadratic 
polynomial in the qi and pi. A system such as this is governed by a symplectic matrix 
M that maps the initial variables Z™ into time-evolved variables Z(t), 

Z{t) = M(t)Z™. (4) 

Let A be given by 

A = hm [MM) , (5) 

where M denotes the matrix transpose of M. The Lyapunov exponents then equal 
the logarithm of the eigenvalues of A [1]. 

It is easy to show that M satisfies the equation of motion (See, for example, Ref. 
[10]) 

where S denotes the symmetric matrix given by 

-1 2m 

H{Z,t) = -Y: S.,Z,Z,, (7) 
«' j'=i 

and where 

J=\ ' ]. (8) 




Here 1 denotes the m X m identity matrix. It follows that the evolution of MM is 
governed by the equation 

4-MM = JSMM - MMSJ. (9) 
dt ^ ' 

Standard methods for obtaining the Lyapunov exponents deal with Af , which is not 
real symmetric (hence the need for reorthogonalization) and which has exponentially 
growing elements. To avoid these difficulties we exploit the fact that M is symplectic 
by making use of the exponential representation of symplectic matrices [6]: Any 
symplectic matrix M can be written in the form 

M = e-'^^e-'^^ (10) 
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where Sa is a symmetric matrix that anticommutes with J and Sc is another sym- 
metric matrix that commutes with J. It is important to note that the second matrix 
on the right hand side of (10) is in fact unitary, so that 

MM = e^'^^\ (11) 

Note that this matrix has fewer degrees of freedom than M and its eigenvectors are 
orthogonaL Rather than attempting to directly integrate Eqn. (9) which would still 
have a large numbers problem, we focus our attention on the exponent JSa in Eqn. 
(10). It is clear that we now do not have a large numbers problem since JSa already 
appears as an exponent. 

To proceed further, one obvious approach is to use an explicit representation 
of exp(JS'a). Such a representation is well known for Sp{2) and has recently been 
found for Sp{4:) (generalizations to Sp{2m) are in progress) [9]. For the purposes 
of this Letter we restrict ourselves to one spatial dimension, where the most general 
symplectic matrix can be written in the form 

— g/i(-B2 cosa+Bs sina)g6_Bi (12) 

where the Bi are basis elements of the Lie algebra sp{2) and where a, h and // are 
real coefficients [6]. It follows that 

jl^jl^ _ g2At(-B2 cosa+Bs sina) {^i) 

Thus, we obtain. 



Finally, it is easily shown that the eigenvalues of this matrix are e^^/*. The Lyapunov 
exponents are then equal to ib/i/t in the limit t — )• oo. With this convenient choice 
of variables, the explicit representation of MM is given by 

(cosh 2 /i + sin a sinh 2 /i cosasinh2/i \ 

. . • (15) 

cosasinh2/i cosh 2/i — sin a sinh 2/i / 

The unknown quantities a and // can grow in time at most as 0{t). We can obtain 
differential equations for these quantities by returning to Fqn. (9), the dynamical 
equation for MM. 

For simplicity, we will assume that H contains no term proportional to qp^ so that 
the matrix in (7) is of the form 

^= I T °V (16) 




After some manipulation, Eqns. (7)-(9) lead to the following: 
dfi 1 

—— = — (S22 — ■Sii I cos a, 

— = sii + S22 — (522 — 511) sinacotli/i. (17) 
at 

From the initial condition Af(0) = /, if we choose /u(0) = 0, then cos^ a(0) = 1, i.e., 
a(0) = or TT. These differential equations form the basis of our method for calculating 
the Lyapunov exponents of Hamiltonian systems: They are stepped forward in time 
numerically till some desired convergence for the exponents, ib/i/t, is achieved. Later 
we will also show how to apply the method to certain non-Hamiltonian systems. 

As a first concrete example, we consider the newly developed "core-halo" model 
which describes beam halo in mismatched charged particle beams [11]. The transverse 
equation of motion for a halo particle in this model, assuming constant external 
focusing, is 

x + x-{l-i]^)f{x,r{t)) = (18) 

where x is the position variable for a halo particle, f[x^r[t)) is the force due to the 
space charge of the beam core, and r[t) is the time dependent rms radius of the core. 
The core radius is assumed to follow the envelope equation 

r + r L__L = o. 19 

Here units have been chosen so that the time independent solution of (19) (i.e., a 
matched beam) is given by r = 1. In these units rj = corresponds to the space charge 
dominated regime, while rj = 1 corresponds to the emittance dominated regime. We 
now assume 

f{^.r) = ^^ (20) 

which has the correct asymptotic behavior: the force is linear when x ^ r, and it 
is inversely proportional to x when x ^ r. The Eqns. (18) and (19) describe a 
driven nonlinear system with a mixed phase space as demonstrated by the strobo- 
scopic plot shown in Fig. 1. The presence of a chaotic band is important because 
particles initially in the core can leak through the broken separatrix and be car- 
ried to large amplitudes. The presence of such large amplitude particles can cause 
unacceptably high radioactivation levels in high intensity linacs planned for future 
accelerator-driven technologies [12]. Leakage through the separatrix can be enhanced 
through particle collisions and recent work has shown that this rate is controlled by 
the Lyapunov exponent [13]. We have computed the Lyapunov exponent for this 

5 



system by integrating (17) with 

S22 = 1 (21) 

where Xq denotes the fiducial trajectory. Fig. 2 displays our result for the Lyapunov 
exponent against time. The slow convergence of the exponent is typical of Hamilto- 
nian systems. 

So far we have dealt with explicitly Hamiltonian systems. However, the only 
real requirement for using our method is that the linearized deviation equations in 
some variables be Hamiltonian. This allows for the inclusion of damped systems in 
our scheme. As an example, we now consider the following general driven nonlinear 
equation 

x + X(l- ex^) X + V'{ x) = a cos[ijjt). (22) 

By appropriate choices of A, e, and V{x)^ this reduces to an assortment of well- 
known equations including van der Pol (A < 0, e = 1, V{x) = (l/2)a;^). Duffing 
(A > 0, e = 0, V{x) = ax^ + /3x'*), and the damped driven pendulum (A > 0, e = 0, 
V{x) = 1 — cos(a;)). In terms of the deviation variable (5, the linearization of (22) 
yields 

^ + A (l - exl) 8 + {V"{xo) - 2eAxoio) = (23) 

where Xq represents the fiducial trajectory. Introducing the new variable A defined 
through 

A = (5e-^(*) (24) 

where 

g = (l - exl) , (25) 

Eqn. (23) reduces to that describing an undamped oscillator with time dependent 
frequency, 

A + (v"{xo) - eXxoio - ^A^ (l - exl)^^ A = 0. (26) 

It is now straightforward to proceed using our method: for linear damping (e = 0) 
the Lyapunov exponents x± of this system are given by 

X± = -^A± limy/io(t) (27) 
where //q follows from solving (17) for the system defined by (26). When e 7^ 0, 

X± = lim y {g{t) ± fioit)) (28) 
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modulo terms that are exponentially suppressed at late times. Figs. 3(a) and 3(b) 
show the Lyapunov exponents of the van der Pol system calculated using our method. 
For the chosen set of parameters our results are in agreement with those of Ref. [4]. 
In contrast with the results shown in Fig. 2, the convergence of the exponents is 
much faster, as is typical of non-weakly damped systems. 

In summary we have described a method for computing Lyapunov exponents that 
exploits the underlying symplectic structure of Hamiltonian dynamics. Just as sym- 
plectic integrators are not a panacea for all time integration problems, we do not 
expect our method to have universal applicability and advantages. However, when 
applicable, the method has certain advantages over standard techniques, most impor- 
tantly the lack of systematic errors associated with intermittent reorthogonalization 
and rescaling. 

The authors acknowledge useful discussions with Alex Dragt, Michael Mattis, and 
Michael Nieto. This work was supported by the U. S. Department of Fnergy at Los 
Alamos National Laboratory and by the Air Force Office of Scientific Research. 
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Figure 1: Stroboscopic plot of the chaotic sea in the core-halo model. Snapshots 
were taken at successive beam minima for 32 test particles. Parameter values were 
r(0) = 0.6, r(0) = 0, and i] = 0.2. 
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Figure 2: Positive Lyapunov exponent for the core-halo model in a typical run. Pa- 
rameters are the same as in Fig. 1. The simulation was run for 10^ periods of the 
driving force. 
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Figure 3: (a) Positive Lyapunov exponent for the van der Pol oscillator with param- 
eters A = —5, a = 5, and u = 2.466 (taken from Ref. [4]). The total time for this run 
corresponds to approximately 10^ periods of the driving force. 
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Figure 3: (b) Negative Lyapunov exponent for the van der Pol oscillator with the 
same parameters as in Fig. 3(a). 
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